Lake Como, nestled in the Lombardy region, is recognized as a crucial freshwater ecosystem esteemed for its biodiversity and ecological significance. Suppose that environmental studies have suggested the presence and spreading of a concerning pollutant within the lake’s waters. The emergence of this pollutant could raise significant concerns about its potential threat to the ecological balance of the aquatic environment and necessitates considerations about its long-term environmental stability. This hypothetical presence prompts a thorough examination to understand its nature, source, and potential implications on Lake Como’s delicate ecosystem.
Assume that several monitoring stations have been strategically established around Lake Como to track and measure pollutant concentrations in its waters. These stations play a pivotal role in continuous surveillance, offering crucial data to assess pollutant levels and their potential impacts on the lake’s ecosystem.
The following interactive figure shows the comprehensive map of Lake Como and includes georeferenced black points denoting the strategic placement of monitoring stations across its expanse. The red points denote potential ‘hotspots’ where, in the hypothetical scenario under investigation, abnormal concentrations of the pollutant have been recorded.
In modeling the spread of the pollutant, it’s crucial to consider the presence of a tributary in the northern region of the lake. Consequently, a transport term has been identified to represent the movement of the pollutant carried by the inflow to other areas within Lake Como. Furthermore, non-homogeneous Dirichlet boundary conditions are applied to boundary regions near the stations that recorded anomalies in the pollutant concentration. Conversely, homogeneous Dirichlet boundary conditions are applied to boundary regions unaffected by the anomalies.
Steady-state solution
In this section, we focus on solving the steady-state problem. Let \(\Omega\) represent our domain of interest. We assume that the concentration of the pollutant, denoted by \(u\), spreads across \(\Omega\) according to the following partial differential equation (PDE):
\[ \begin{cases} -\mu \Delta u + \boldsymbol{\beta} \cdot \nabla u = 0 \qquad & in \ \Omega \\ u = g|_{\partial \Omega} \qquad & on \ \partial \Omega, \end{cases} \] where \(\mu \in \mathbb{R}\), \(\boldsymbol{\beta} \in \mathbb{R}^2\) are the diffusion coefficient and the transport coefficient, respectively. The function \(g\) used to enforce the boundary conditions consists of the sum of two Gaussian-like functions centered at the locations of the stations that detected anomalous concentrations of the pollutant.
The upcoming sections will illustrate how to obtain a discrete
solution of the problem using femR package.
library(sf, quietly = T)
library(mapview, quietly = T)
library(femR, quietly = T)
library(rmapshaper, quietly = T)
lake_bd <- read_sf("data/deims_sites_boundariesPolygon.shp")
lake_bd_simp <- rmapshaper::ms_simplify(lake_bd, keep= 0.5,
keep_shapes=T)
st_crs(lake_bd_simp) <- 4326The femR package includes utilities to read a
sfc (Simple Features Collection) object, such as the
boundary of the Lake, and generating meshes by specifying the
maximum_area of each triangle and the
minimum_angle between adjacent sides of the triangles.
# create the physical domain
domain <- Domain(st_geometry(lake_bd_simp))
# create the mesh
mesh <- build_mesh(domain, maximum_area = 0.25e-5, minimum_angle = 20)Furthermore, the Mesh class overloads the
st_as_sfc method, enabling visualization of the mesh object
on an interactive map using the mapview package, as
demonstrated in the following code snippet.
mesh_sf <- st_as_sfc(mesh)
map.type <- "Esri.WorldTopoMap"
mapview(lake_bd_simp, col.regions = "white", lwd = 2,
alpha.regions=0.2, legend=F, map.type=map.type) +
mapview(mesh_sf, col.region="white", alpha.regions=0.2, legend=F, alpha=0.4,
lwd=0.4, map.type=map.type) The next step involves defining the FunctionSpace where
to find the solution of the PDE. By default, first-order finite elements
are considered during the construction of Vh. Then, the
solution of the problem is instantiated as an element of
Vh.
Vh <- FunctionSpace(mesh)
u <- Function(Vh)Then, we define the differential operator in its strong formulation as follows:
Both the forcing term of the differential problem and the boundary
condition can be defined either as standard R functions or
as simple numerical values, particularly when representing constant
values.
# Dirichlet boundary conditions
g <- function(points){
res <- matrix(0, nrow=nrow(points), ncol=1)
pts_list <- lapply(as.list(as.data.frame(t(points))), st_point)
dists2 <- st_distance(st_geometry(sources),
st_sfc(pts_list, crs=st_crs(sources)))
dists2 <- matrix(as.numeric(dists2),
nrow=length(st_geometry(sources)), ncol=nrow(points))/(2.5*10^3)
for(i in 1:length(sources))
res <- res + Q*exp(-dists2[i,])
return(res)
}Finally, we build a Pde object passing the differential
operator L, as first parameter, and the forcing term, which
is a constant in this specific case, as second parameter:
# build PDE object
pde <- Pde(Lu, 0.)
# set boundary conditions
pde$set_boundary_condition(fun= g, type= "dirichlet")We can compute the discrete solution of the problem calling the
solve method:
pde$solve()Finally, it is possible to visualize the computed solution on an
interactive map using the mapview package as the following
code snippet shows.
coeff <- apply(mesh$get_elements(), MARGIN=1,
FUN= function(edge){
mean(u$coeff[edge])
})
U <- st_sf(data.frame(coeff = coeff), geometry= mesh_sf)
mapview(U, alpha.regions=0.8, alpha = 0, map.type=map.type)Time-dependent solution
Now, we aim to estimate the evolution of pollutant concentrations in Lake Como. We consider a Gaussian-like function \(g\) from the previous section, appropriately modified to account for the presence of time. This function will be used to define both the initial condition for the pollutant concentration and to impose the boundary conditions. Therefore, the partial differential equation modeling the evolution of the pollutant concentration in Lake Como is as follows:
\[ \begin{cases} \frac{\partial u}{\partial t} -\mu \Delta u + \boldsymbol{\beta} \cdot \nabla u = 0 \qquad & in \ \Omega \times (0,T) \\ u = g & in \ \Omega, \ t=0 \\ u = g|_{\partial \Omega} \qquad & on \ \partial \Omega \times (0,T), \end{cases} \]
The following window wraps all the steps needed to estimate the
evolution of the pollutant concentration relying on femR
package.
time_interval <- c(0., 0.3)
times <- seq(time_interval[1], time_interval[2], length.out=50)
# function space
Vh <- FunctionSpace(mesh %X% times)
# PDE solution
u <- Function(Vh)
# differential operator
Lu <- dt(u) - mu*laplace(u) + dot(beta, grad(u))
# dirichlet bc
g <- function(points, times){
res <- matrix(0, nrow=nrow(points), ncol=length(times))
pts_list <- lapply(as.list(as.data.frame(t(points))), st_point)
dists2 <- st_distance(st_geometry(sources),
st_sfc(pts_list, crs=st_crs(sources)))
dists2 <- matrix(as.numeric(dists2),
nrow=length(st_geometry(sources)), ncol=nrow(points))/(2.5*10^3)
for(i in 1:length(sources)){
res[,1:length(times)] <- res[,1:length(times)] + Q*exp(-dists2[i,])
}
return(res)
}
# initial condition
u0 <- g(mesh$get_nodes(), times[1])
# build PDE object
pde <- Pde(Lu, 0.)
# set initial conditions
pde$set_initial_condition(u0)
# set boundary conditions
pde$set_boundary_condition(fun= g, type= "dirichlet")
# compute the discrete solution
pde$solve()
# plot
plot(u)